Thermal fluctuations in the lattice Boltzmann method for non-ideal fluids 
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We introduce thermal fluctuations in the lattice Boltzmann method for non-ideal fluids. A 
fluctuation-dissipation theorem is derived within the Langevin framework and applied to a spe- 
cific lattice Boltzmann model that approximates the linearized fluctuating Navier-Stokes equations 
for fluids based on square-gradient free energy functionals. The obtained thermal noise is shown 
to ensure equilibration of all degrees of freedom in a simulation to high accuracy. Furthermore, we 
demonstrate that satisfactory results for most practical applications of fluctuating hydrodynamics 
can already be achieved using thermal noise derived in the long wavelength-limit. 



I. INTRODUCTION 



Owing to its flexibility and easily parallelizable nature, 
the lattice Boltzmann (LB) method has by now become 
an established tool for solving the Navier-Stokes equa- 
tions for simple as well as complex fluids The sim- 
ulation of systems with phase coexistence is not only 
important for applications, such as wetting and thin- 
films 0, 01, but also interesting from a theoretical point 
of view [J, It is well known that the thermal mo- 
tion of the fluid particles becomes relevant already be- 
low the micro-scale, leading — for example — to the Brow- 
nian motion of suspended solid particles 6 8]. But also 
in pure fluid systems, thermal noise plays an important 
role close to phase transitions @ or hydrodynamic in- 
stabilities EH , and has recently been shown to also 
have significant effects on nanoscopic free-surface flows 
Simulation of such behavior with a determinis- 
tic method, such as LB, requires the inclusion of explicit 
noise sources in the underlying equations. In this work, 
we will discuss how thermal noise can be modeled within 
the LB method for non-ideal fluids. 

The history of thermal fluctuations in the Boltzmann 
equation dates back to Kadomtscv (l5| . who first ap- 
plied the Langevin approach to the Boltzmann equa- 
tion of a dilute gas. It was shown later by Bixon and 
Zwanzig [l6| . and independently by Fox and Uhlenbeck 
[l7j|, that this approach in fact leads to the well-known 
equations of fluctuating hydrodynamics fl8j in the limit 
of large length and time scales. Generalizations of the 
Boltzmann-Langcvin equation to non-ideal gases have 
been discussed by Klimontovich It has been shown 
by Kim and Mazenko [2(|, that the expressions for the 
fluctuating stress tensor known for a simple fluid essen- 
tially remain valid also for a fluid described by a square- 
gradient free energy functional. 
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In the context of the LB method, so far, only the 
fluctuating ideal gas model has been studied systemat- 
ically, starting with the work of Ladd (U [22j ■ There, 
the Landau-Lifshitz theory of fluctuating hydrodynam- 
ics [lH was implemented by adding a fluctuating compo- 
nent to the LB stress modes. While the model satisfied 
the fluctuation-dissipation theorem (FDT) at the hydro- 
dynamic level, it was soon realized that it failed to give 
full equilibration of momentum [23j . As first pointed 
out by Adhikari et al. J24[, in order to ensure correct 
equipartition of fluctuation energy at all length scales, 
noise must not only be added to the stresses, but to all 
dissipative modes that exist for a given model. This was 
confirmed subsequently by Diinweg et al. [25| using a lat- 
tice gas analogy. Notably, Dufty and Ernst [26| gave the 
first general treatment of LB-Langcvin models, including 
a derivation of the appropriate FDT. We finally mention 
that there also exist a number of finite-volume schemes 
for the direct integration of the fluctuating Navier-Stokes 
equations [27j . 

In the present work, we study thermal fluctuations in 
a non-ideal fluid within the Langevin framework. First, 
the fluctuating hydrodynamic equations for a fluid based 
on a square-gradient free energy functional and the asso- 
ciated FDT are discussed for a continuum system. Next, 
we present a derivation of the FDT appropriate to the 
non-ideal fluid LB method. Consistency requires that the 
stochastic LB equation (LBE) leads to the same form of 
the fluctuating stress tensor that is required by the FDT 
derived independently at Navier-Stokes level. The theory 
is then applied to the modified-equilibrium model of Swift 
et al. [28|, [29[. We find that, for this model, the noise 
must in general be spatially correlated to ensure ther- 
malization at all length scales. Additionally, we demon- 
strate that, at least for certain regions in the parameter 
space, satisfactory results at large length scales can also 
be achieved using an approximate, spatially uncorrelated 
form of noise derived in the hydrodynamic limit. Finally, 
it is shown that capillary fluctuations can successfully 
be simulated using uncorrelated noise in the modified- 
equilibrium model. 
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FIG. 1: (Color online) Thermodynamic model: (a) bulk pres- 
sure po and shape of the free energy density fo (not to scale) , 
(b) sound speed squared. Parameters: pv = 0.1, pL = 1-0, 
ft = 0.01. 



II. CONTINUUM THEORY 

We first review the general physical background of the 
non-ideal fluid models considered here. In this work, the 
following convention for the spatial and temporal Fourier 
transform of a quantity a(r,t) is applied: a(k, u)) = 
(l/27r) d / 2 / drdta(r,t)exp(i(k-r-ut)), where d is the 
spatial dimension. 



A. Thermodynamics 

Our treatment is based on non-ideal fluid models de- 
scribed by a square-gradient free energy functional, 



Hp] 



dv(fo(p) + ^ P f) • (1) 



Here, fo is the bulk free energy density and k is the 
'square-gradient' parameter, which can be related to the 
surface tension and interface width. For fo, we take a 
simple Landau-type double- well potential [H, HO, HH : 



f (p) =P{p- pv?{p- PL? 



(2) 



where pv,L are the desired equilibrium vapor and liquid 
densities and the parameter f3 is inversely proportional 
to the compressibility of the fluid. The bulk pressure po 
(equation of state) can be computed from the free energy 
density via po = pd p fo — fo- Fig. [1] illustrates the typical 
shape of the bulk free energy, bulk pressure and speed of 
sound, c 2 s = dp /dp = pd 2 f /dp 2 . 

In a single phase, the parameter (3 is related to the 
speed of sound c s by 



where po corresponds to cither pl or py, depending on 
which phase c s is referring to. Typical values of c s in our 
simulations range from 0.04 to 0.3 in lattice units (l.u.) 
(see section W\ . hence the compressibility of the simu- 
lated non-ideal fluid is strongly enhanced compared to 
the ideal gas case, where c Si idcai = \/l/3 ~ 0.57 (Fig. [1]). 
Coexisting phases in a square-gradient fluid are generally 
separated by a diffuse interface [10, H|[ , which — for the 
above form of the free energy potential — has a width of 



1 



P PL ~ PV 

The surface tension follows as 

(PL - Pv) 3 



(3) 



(4) 



The thermodynamic pressure tensor that follows from 
the above free energy functional is given by j32hj36l | 



P 



ih 



Po 



npV z p--\Vp\ z )l + K^Jp)®{Vp). (5) 



Due to the square-gradient term, the pressure tensor re- 
ceives non-local contributions in addition to the bulk 
pressure po- In the Navier-Stokes equations, the diver- 
gence of this tensor appears, which can be written as the 
sum of a gradient of the thermodynamic bulk pressure 
and a force-like term: 



V • P th = Vp - np\JV z p. 



(6) 



Another route to derive the thermodynamic interac- 
tion force consists of directly computing the effective 
chemical potential from the free energy functional, 

-2 

(J, = — = p a - K V p , 
dp 

with po = dpfo- The effective "chemical" body force that 
is acting on a fluid element is thus given by 



F cff = -pV^t = -Vpo + tpVV 2 /9 = -V • Pi 



2>Po(PL ~ Pvf 



B. Fluctuations 

We consider fluctuations around a quiescent, homo- 
geneous equilibrium state of density p and vanishing 
macroscopic flow velocity u = 0, i.e. p(r) = p + Sp(r) 
and u(r) = <5u(r). The fluid momentum is given by 
<5j = poSu. 

Density fluctuations. The essential difference between 
an ideal gas and a non-ideal fluid is the fact that, in the 
latter, density fluctuations are spatially correlated. The 
density correlation function (static structure factor) can 
be determined by expanding the free energy functional up 
to second-order in the density around equilibrium [0, H3, 
|38| . By this, one obtains a Gaussian probability density 
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in Fourier-space, with a variance given by an Ornstcin- 
Zcrnike type structure factor 



(<Jp(k)*p(k')> = S(k)<f(k + k') = 



p k B T 



p K 



k 2 



S(k + k') 



(7) 

Here, k B the Boltzmann constant and T is the tempera- 
ture of the fluid. The correlation length associated with 
the density fluctuations is given by ^/pqk/c s , which is 
directly proportional to the interface width, cq. ([3]). 

Momentum fluctuations. In an equilibrium fluid, the 
momenta of the fluid particles are always uncorrelated 
[37j . hence, by equipartition, the equal-time momentum 
correlation function is given by [l8j 



{6j a (k)5jp(k')) = Po k B T6 aP S(k + k') 



C. Hydrodynamics 



(8) 



Using expression © for the divergence of the pressure 
tensor for a fluid with an underlying square-gradient free 
energy functional, one obtains the linearized stochastic 
Navier-Stokes equations in d dimensions as 



d t Sp = -poV ■ u, 
p d t u = -c 2 VSp + np V(V 2 5p) 

+ nV 2 u + (( + r)[l - 2/d]) V(V ■ u) - V • R. 



(9) 



Here, 77 and £ are the shear and bulk viscosities and R is 
the random stress tensor [l8| • We now Fourier transform 
the space dependence of dp, u and R, and separate the 
velocity and the random stress tensor into longitudinal 
and transverse components, u = u/k + u t , where k = 
k/|k|, ui = u • k, and Ut = u • (I— kk), and analogously 
k • R = Rik + R t , where R ( = k ■ R ■ k, and R t = 
k ■ R ■ (I - kk). Wc thus arrive at [H M, 53] 



d t Sp = ipokui 



(10) 



dm = ik (c 2 + p Kk 2 ) — - vik 2 ui + —Ri , (11) 

Po 



A) 



d t u t = -v t k u t H R t . 

Pa 



(12) 



where we have introduced the longitudinal and transverse 
kinematic viscosities vi = [£ + n (2 — 2/d)]/ po, and v t = 
V/Po- 

The essential observation is that the above equations 
are identical to those for a simple bulk fluid (i.e., a fluid 
where k = 0), if one introduces a wavenumber dependent 
sound speed 

c 2 (k) = c 2 + p oK k 2 = p k B T/S(k) (13) 

in the latter. Further, we note that also for a square- 
gradient fluid it remains true that density fluctua- 
tions only couple to longitudinal momentum fluctuations, 



while transverse momentum fluctuations are completely 
decoupled from the other variables. 

The fluctuation-dissipation theorem relates the equal- 
time correlation function of the fluid momentum to 
the correlation function of the fluctuating stress ten- 
sor R, which is assumed to obey a Gaussian prob- 
ability distribution and have a Markovian character, 
{R a p{ri,ti)R a p{r2,t 2 )) = A a p(r-y - r 2 )S(t 1 - t 2 ), with a 
variance A a p to be specified below. Since the non-ideal 
fluid interactions enter the hydrodynamic equations only 
through a modified speed of sound, one might expect 
that the expressions of the fluctuating stress tensor of a 
square-gradient fluid and a simple bulk fluid are identi- 
cal. This is in fact true [23], and it holds even in the case 
of non-linear fluctuating hydrodynamics, if one properly 
takes into account the local values of the transport coeffi- 
cients in the random stress tensor (which then represents 
multiplicative noise) . 

In the following, we explicitly demonstrate this fact 
and derive the expression for the variance of R within the 
Langcvin framework. Inserting the continuity equation 
(fTU|) into cq. (|TTj) and solving for the longitudinal velocity, 
one obtains 

ik 

d 2 ui = -k 2 c 2 s (k)ui - uik 2 d t ui H d t Ri ■ 

Po 

Fourier transforming in time, we obtain the linear re- 
sponse relation: 



po (oj 2 — k 2 c 2 (k) — iujvik 2 ) 

= Xi(k,uj)Ri(k,oj), (14) 



keeping in mind that the longitudinal susceptibility xi ( w ) 
has poles at complex frequencies. Squaring this equation, 
averaging over the noise and employing the white-noise 



property of R, (l-Ro^k, £ 



^4 Q( a(k), the equal-time 



correlation function of it; follows after an inverse Fourier 
transform as 



MM = o)| 2 ) = i- J du{\ Ul {k^ 



2tt 



duj\ X i(k,oj)\\ (15) 



Here, we introduced the longitudinal component Ai of the 
variance A a p. The integral over Xi( w ) can De computed 
using contour integration, giving J du\xi{uj)\ 2 = 7r/ p 2 vi- 
In order to obtain a thermally equilibrated fluid, the noise 
average of the velocity correlator in (fl~5j) is set equal 
to the thermal average. Since equipartition demands 
(|M/(k)| 2 ) = k B T/p , we finally obtain 



(\R[(k,uj)\ 2 )=2k B Tpovi. 



(16) 



The corresponding relation for each transverse compo- 
nent of R follows analogously as 



(|i?t(k,w)| 2 ) =2k B Tp Q v t . 



(17) 
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Hence, the "classical" FDT of fluctuating hydrodynamics 
for a simple bulk fluid is recovered, with a random stress 
tensor that is uncorrclated in space and time. Fourier- 



J 



(R a p{T,t)R yS (r',t'))=2k B T 



transforming back to the real space and time domain, the 
FDT assumes the well-known form [H, Eol 



5{r - r')S(t - t') 



(18) 



Note that the last expression is valid only if the viscosi- 
ties are independent of k, as is the case for a simple bulk 
fluid. The particular non-ideal fluid LB model we con- 
sider in section lTlI C[ however, entails a fc-dependent bulk 
viscosity [see eq. (|35|) ] and hence the random stress tensor 
becomes spatially correlated in this case. 

In an inhomogeneous state, i.e., in the presence of in- 
terfaces, the fluid is governed by non-linear equations of 
motion since the non-linear terms originating from the 
pressure tensor ([5]) can not be neglected anymore. How- 
ever, small fluctuations around a solution of the full non- 
linear equations can still be locally described by the same 
linearized equations of motion ([§]) as in the uniform case. 
Therefore, the FDT of fluctuating hydrodynamics (fT8|) 
is expected to remain valid if the effect of the inhomo- 
gencity is taken into account in the local values of the 
thermodynamic quantities and the transport coefficients 
2C , . This renders the noise effectively multiplicative 
42(| . Applying the linear Langevin formalism to general 
non-equilibrium situations can often be justified along 
similar arguments after a local-equilibrium assumption 
has been made I4fl, [431 . 



III. LATTICE BOLTZMANN MODELING 

We now turn to the LB model that corresponds to the 
continuum theory of non-ideal fluids presented in the pre- 
ceding section. Although we focus on a two-dimensional 
system, all derivations arc kept as general as possible 
and are easily applied to three dimensions. The fluc- 
tuating LBE is most conveniently developed in terms of 
moments of the distribution function. The moment space 
is divided into the conserved, transport and ghost (or ki- 
netic) sectors. The FDT is derived below in Fourier-space 
by treating all LB modes on an equal footing, which, as 
first pointed out in 24j, is a necessary prerequisite to 
achieve complete thermalization of the fluid. The eval- 



uation of the FDT requires, besides information on the 
relaxation and interaction behavior (which is provided by 
the LB model itself), additional information in the form 
of the correlation matrix of the LB modes. The latter in- 
gredient has to be derived from a statistical mechanical 
framework. The FDT is finally applied to the modified- 
cquilibrium non-ideal fluid model of Swift et al. [28|, H^] . 
Notation and conventions. The spatial Fourier trans- 



form of a quantity a(r) defined on the lattice is computed 
according to a(k) 



l 



J2 r e lk r a(r) , where n is the total 
number of lattice points in the system. Since a is usually 
a real quantity, it is sufficient to consider just the first 
quadrant of the first Brillouin zone, k a = . . .it, where 
k a = tt corresponds to a physical length of 2 l.u. Fourier 
transforming derivative operators on a lattice requires 
additional care, as described in appendix |Bj In general, 
Greek indices refer to Cartesian coordinates, while Latin 
indices refer to the lattice directions or the LB moments. 
Repeated free indices are to be summed over. In the fol- 
lowing, the quantity a s = y/T/3 is a constant specific to 
the chosen lattice and agrees with the speed of sound of 
the ideal LB gas. It has to be distinguished from the 
actual speed of sound c s of a non-ideal fluid. 



A. Introduction 

We begin by reviewing the necessary theory of the de- 
terministic LBE, which is given by 



fi(r + c h t+l) = fi+Aijifj 



IP) 



■Ft 



(19) 



Note that the r- and ^-dependences have been suppressed 
on the right hand side of the equation. Here, Fi describes 
a possible body force and Ay is a general matrix relax- 
ation operator [3, |4j| . The equilibrium distribution f^ q 
is model-dependent and will be specified later. In this 
work, we consider a D2Q9 lattice, hence i = 1,...,9. 
For a general discussion of the LB method for simple 
and complex fluids, we refer to the literature (ll. 14114491 ] . 

For the present purposes, it proves to be most conve- 
nient to work in the space of moments m a (a = 1, . . . , 9) 
of the distribution function fi (Hoj . This can be achieved 
by constructing a set of orthogonal basis vectors T a i from 
the lattice velocities Cj. Orthogonality is measured with 
respect to the weighted scalar product [24| . 



(Ta\T b , 



WiT r ,,T, 



hi 



Na&ab 



(20) 



where N a is the length of the ath basis vector T a , N a — 
Si w i^ai- The weights Wi are identical to the ones used 
in the definition of the (ideal gas) equilibrium distribu- 
tion. For a standard D2Q9 lattice, these are 



Wi = 4/9, W2...5 = 1/9, w 6 . 



1/36. 
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With this choice, the projection of the (ideal gas) equi- 
librium distribution onto the ghost modes is eliminated 
[2~H |. Although this is in principle not necessary for the 
present developments, we will adopt this choice hence- 
forth, since it emphasizes the physical background of the 
model and is computationally advantageous. 

The moments are now defined as the projections of the 
distribution function onto the basis vectors, 

In turn, any vector defined in velocity space (such as the 
distribution function fi) can be expanded in terms of the 
orthogonal basis vectors, 

fi = (T~ 1 )i a m a = w l T ai m a /N a . 

In this relation, the expression for the inverse transfor- 
mation matrix, (T -1 )^ = WiT a i/N a , has been used. 

Once a suitable basis set is chosen (see below) , one con- 
structs a collision operator A that is diagonal in moment 
space by setting A = T _1 AT, where A = diag(A a =i,...,g) 
is a diagonal matrix of relaxation parameters A a . Hence, 
the basis vectors T a are eigenvectors of the generalized 
collision operator A with eigenvalues A a . The A a can 
be expressed in terms of relaxation times r a through 
A a = — 1/t . Well-known stability requirements impose 
the restriction r a > 1/2 [lj. Rewriting the right hand 
side of eq. p^|) in terms of moments, the LBE becomes 

fi(x + a,t+l) = T- 1 [m a + X a {m a - m^) + mf ] , 

(21) 

where mf = T a iFi are the moments of the forcing term. 

The D2Q9 basis set used in the present work is sum- 
marized in Table Q] [25j . The first three rows cover the 
conserved hydrodynamic moments, i.e. the density and 
momentum 



The next three rows contain the non-conserved hydrody- 
namic (transport) moments, and the last three contain 
the ghost (or kinetic) moments. The moment e describes 
a pressure or bulk stress mode, with an eigenvalue Af, re- 
lated to the bulk viscosity. p ww and p xy are shear modes, 
with a common eigenvalue A s related to the shear viscos- 
ity. The ghost modes consist of a ghost density mode 
e and ghost vector current q a , in line with the duality 
prescription (5lj . Using the numerical expressions for 
the lattice velocities, the transformation matrix for the 
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N a 
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1/3 


3x 
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Ciy 


1/3 
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4/9 


PwW 
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Pxy 
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(3c? - 4)c„ 
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Qx 
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(3c? - 4)c lH 
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9cf - 15c? + 2 


16 
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A £ 



TABLE I: Basis set of the D2Q9 model. T ai denotes the basis 
vector, N a its length, m a is the designation of the correspond- 
ing moment and A a is its eigenvalue in the collision operator. 
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B. FDT 

The fluctuating LBE is obtained from the determinis- 
tic LBE, cq. ((2~T|) . by adding random noise variables £ a to 
the collision step. A small fluctuation around a uniform, 
global equilibrium state of density po an d vanishing flow 
velocity, Sfi(r,t) = fi(r,t) - f^(po,u = 0), or equiva- 
lently, <5m (r, t) = m a (r, t) — m® q (/9o, u = 0), then evolves 
according to the linearized equation 

8fi(T + c h t + l) = Tr 1 ^ + \ a )Sm a - X a 8m^ 

+ <5mf + £ Q ]. (23) 

The noise is assumed to be uncorrected in time and 
drawn from a Gaussian probability distribution, whose 
covariance matrix can be expressed as (£ n (r, t)^b(r', t')) = 
S a fa(r — r')St,r , assuming translational invariance. It fur- 
ther is assumed that the £ a are uncorrclated with the LB 
modes Sm a . The matrix S a b will be determined below 
by means of the fluctuation-dissipation theorem. 

In order to rewrite the fluctuating LBE (|23|) fully in 
terms of moments, we apply a spatial Fourier transform 
and introduce the Fourier-transformed advection opera- 
tor in moment space by [52j 

A ab (k) - T aj exp(-ik • c,)(T" 1 ), b . (24) 
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This operator has the important property that A~ h x (k) = 
A* b (k) , where * denotes complex conjugation. Collecting 
the effect of relaxation and interactions in a linearized 
collision operator O(k), we finally obtain 

Sm a (k, t + l) = A ab (-k) [(I + n) bc Sm c {k, t) + &(k, t)} . 

(25) 

The computation of f2 requires only the knowledge of 
(5m oq and 8m F and will be performed in section fill CI for 
a specific LB model. 

The derivation of the FDT for a non-ideal fluid model 
proceeds analogously to [13, Hlj]. We multiply eq. (|25|) 
with 8md{— k, t + 1) from the right, average over the 
noise distribution and assume stationarity of equal time 
correlators. Introducing the equal-time correlation ma- 
trix of the modes as G a b(k) = (5m a {k)5mb(— k)) and the 
(Fourier-transformed) covariance matrix of the noise as 
S ab (k) = (£ a (k)£f,(-k)), we obtain the FDT in the form 



3(k) 



A(k)G(k)A(-k) T - [l+fi(k)] G(k) [l+0(-k)l T . 

(26) 

Note that, due to the absence of non-linearities, the above 
expression can be independently evaluated for each point 
in fc-space. 

The equilibrium correlations of the modes, represented 
by the matrix G, are not immediately provided by the LB 
scheme, but must instead, as for any Langevin equation, 
be derived with the help of a statistical mechanical theory 
p3 - l26^ . A suitable ansatz for a general non-ideal fluid 
is provided by the following relation for the equal-time 
correlations of fluctuations in the distribution function: 

{SfiMSfjQi!)) = [/i/i[5(k)/p - A i]/p + /i/ i <5 ij ]5 kl _ k /. 

(27) 

Here, fi = f^ q (po,u = 0) denotes the global Maxwellian 
of the uniform reference state and 5(k) is the structure 
factor of the non-ideal fluid. The parameter /i can be 
interpreted as the mass of a fictitious fluid particle and 
must be determined such that eq. (|27p leads to the correct 
expression ([8]) of the equal-time momentum correlations. 
The correlation matrix G is obtained from ([27| through 
a basis transformation, 



G ab (k) = T a JT bj (6f i (k)6f j (-k)) 
= fh a fh b [S(k)/p - /x]/p 



(28) 



where fh a = T a ifi. A derivation of relation (|27p from 
continuum kinetic theory is presented in appendix [X] 
The above relation is expected to be appropriate to 
non-ideal fluid models that are based on the ideal gas 
equilibrium distribution, which is the LB equivalent of 
the Maxwellian distribution used in continuum kinetic 
theory. Conversely, non-ideal fluid models employing 
a non-standard equilibrium can be expected to require 
a different ansatz from (|2"7|) . This is indeed the case 
for the modified-equilibrium model considered below [see 
eq. ([52")) ]. The general theoretical status of (|2"7| in the 
context of the LB method for non-ideal fluids will be fur- 
ther investigated in future works. 



In the hydrodynamic regime, i.e. at small k, the vari- 
ances of the noise variables pertaining to the stress modes 
(here £4,5,6) can be determined directly using the FDT 
of fluctuating hydrodynamics (appendix [C]). This fact 
allows for an independent check of the noise constructed 
on LB level, eq. (|26]l. In many cases it turns out that 
satisfactory simulation results can already be achieved 
by using noise evaluated for k — > 0, which is then spa- 
tially uncorrelated by construction, similarly to the ideal 
gas case. This behavior can be expected to apply when- 
ever the non-ideal fluid interactions enter the dynamic 
equations in a fully reversible way — and hence, not give 
additional contributions to the dissipative terms — or if 
the irreversible contribution is sufficiently weak. Since 
A = I at k = 0, the noise covariance matrix reduces in 
the zero wavelength-limit to 



3(0) 



lim 3(k) 



GVl 1 -nG-VLGtt 1 



(29) 



where all quantities on the right hand side are evalu- 
ated for k — > 0. This uncorrelated form of the noise can 
be constructed independently on each lattice site in real 
space. 

Moreover, as remarked in section Hi C| one can usually 
assume the FDT of fluctuating hydrodynamics (|T5|) . de- 
rived from the linearized Langevin equations, to remain 
still valid even for an inhomogeneous fluid, if the local 
values of the thermodynamic quantities and the trans- 
port coefficients are taken into account in the computa- 
tion of the noise. This implies, that uncorrelated noise 
described by 5(0) can be readily applied to inhomoge- 
neous systems, as will be further explained in the next 
section. 



C. Application to a non-ideal fluid model 

We now apply the general FDT derived above to the 
modified equilibrium model proposed by Swift et al. [HI, 
I29I |53| . In this approach, the non-ideal fluid interactions 
are derived from a square- gradient free energy functional, 
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- Kp v 2 p) - 3« \(d x p) 2 + (d y p) 2 ] 



TABLE II: Equilibrium moments of the modified-equilibrium 
non-ideal fluid LB model. The C n denote Galilean-invariance 
correction terms. 
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cq. (HJ), and hence, the stochastic version of this model is 
supposed to approximate the fluctuating Navier-Stokcs 
equations (0 on large length and time scales. The non- 
ideal fluid interactions enter this model through a stress 
contribution to a modified equilibrium distribution whose 
moments are stated in Table [TTJ The C n are correction 
terms that ensure Galilean-invariance of the model up 
to 0(Ma 2 ) [54| . Since they are at least of second order 
in p und u, their specific form is not important for the 
linearized model. We notice, that, in contrast to the 
ideal gas model, the equilibrium distribution has a non- 
vanishing projection onto a ghost mode. The evolution 
equation for the modified-equilibrium model is given by 
eq. (|2"Tj) . without a forcing term (i.e. = 0). 

In the linearized model, after applying a spatial Fourier 
transform (see appendix iBl) and writing Spo = c 2 s 5p, the 
fluctuations of the equilibrium modes follow from Table [TT1 
as 



where the dots indicate zeros for short. 



To compute the equilibrium correlation matrix G = 
{Sm a Sml} = T a iTi,j{5fiSf*} for the modified-equilibrium 
model, we propose here the ansatz 



(6 f t (\i)5f ,(-*)) 



S(k) 
Po 



fi(k)6 i: 



(32) 



5rrC = {Sp, Sj x ,Sj y ,d(k)Sp, 0, 0, 0, 0, -d(k)Sp} , (30) 

where we invoked the definition of the generalized speed 
of sound, c 2 (k) = c? s + ponk 2 , and defined d(k) = 
6 [c 2 (k) — erf] . The model can eventually be brought into 
the general form of eq. (|25p by introducing the matrix 
collision operator as 





f " 




. . . \ 




fi(k) = 


-A 6 d(k) . . 


X b ■ ■ 

. X s . 

■ ■ x s 




, (31) 

I 




{ A e d(k) . . 




x q . . 
■ x q . 
. ■ A J 





where S(k) is the structure-factor, defined by eq. (J7J, 
taking into account the proper lattice derivative opera- 
tors (see appendix iBj) . and /i(k) = fP(k 7 u = 0) is the 
equilibrium distribution of the model evaluated for van- 
ishing flow velocity. This expression can be interpreted as 
the natural generalization of the corresponding ideal gas 
relation, (SfiSfj) = Sidfi$ij/po, to the present non-ideal 
fluid model. The connection of relation (|32[) to the con- 
tinuum kinetic theory of non-ideal fluids and further mo- 
tivations pointing to its validity are detailed in appendix 
El In moment space, the correlation matrix becomes 



G(k) 



S(k) . . 

fa 2 . 
• fa 2 


5(k)d(k) 


-S(k)d(k) \ 


5(k)d(k) 


N 4 S(k) 


25(k)d(k) 




N 5 f . 






• N 6 f 




-S(k)d(k) . 


25*(k)d(k) . 


N 7 f . 

■ N s f 

. N 9 (S(k)+3T)/4j 



(33) 



where the N a are the lengths of the basis vectors and T = poksT/a 2 for short. 

Evaluating the noise covariance matrix (j2"o| using the expressions for f2(k), eq. (|3T1) . and G*(k), eq. (|3"3"]l . shows that 
the advective contribution cancels, i.e. AGA^ = G for all k, just as in the ideal gas model [24j and in the continuum 
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Boltzmann-cquation [19( . The noise covariance matrix is finally obtained as 





( ■ 




\ 

■ 


H(k)= pok f 




Ni [2 - 3c 2 (k)] A e 

N 5 \ s . 
■ N 6 X S 


. 12[ C 2(k)-(T 2 ] A e£ 




{. . . 


12 [^(k)-^ 2 ] A ee . 


N 7 X q . 
■ N 8 X q 

. iV 9 [|- fc 2 (k)] A e / 



where we defined A a = A a (2 + A a ) and A ee = A e + A e + 
A e A e for short. This result shows that the modificd- 
cquilibrium model requires spatially correlated noise to 
satisfy the FDT at all scales. This is in complete agree- 
ment with the FDT of fluctuating hydrodynamics (see 
appendix [C]) for this model: the bulk viscosity of the 
modificd-cquilibrium model is found to be given by 



C(k) = p a 2 n 



;(k) 



(35) 



and hence (in contrast to an ideal gas model) is 
wavenumber-dependent [46[ . Indeed, we recognize in S44 
the presence of the same correction factor 2 — 3c 2 (k) that 
also appears in the bulk viscosity. The shear viscosity is 
the same as in an ideal gas model, r\ — p^a 2 (r s — i), in 
further agreement with result (|34|) . We shall take these 
observations as crucial hints to the correctness of the 
ansatz (|3"2"|) for the modificd-cquilibrium model. Results 
([3~4"|) and ([33)1 indicate that the non-ideal interactions are 
implemented in this model in a way that is not fully re- 
versible. 

A useful approximation to the full noise matrix consists 
of evaluating S in the limit k — > 0, resulting in spatially 
uncorrelated noise. Although this form of noise satis- 
fies the FDT of fluctuating hydrodynamics strictly only 
in the zero wavenumber-limit, simulations indicate that, 
at least for certain parameter ranges, satisfactory ther- 
malization is obtainable also well in the finite k regime. 
This can be explained by the weak /c 2 -dependence of 
the generalized speed of sound c s (k), and hence of the 
noise (|3~4")l and the bulk viscosity (|3"5"|) . Noise defined by 
S(0) = lini/^o 2(k) differs from the noise of an ideal gas 
model by the presence of cross-correlations between stress 
(e) and ghost (e) noises, which originate from the cross- 
correlations in the equilibrium correlation matrix Q33p . 
For the simulation of hydrodynamics at small wavenum- 
bers, these cross-correlations are immaterial as the ghost 
modes are decoupled from the hydrodynamic modes. 

The requirement of positive-definiteness of H imposes 
a restriction on the allowed values of the square-gradient 
parameter k and the sound speed c s through their rela- 
tion to the generalized speed of sound c 2 (k) = c 2 + poKk 2 



[65l| . This is demonstrated in Fig. [31 where the white 
region corresponds to parameter combinations that en- 
sure a positive-semidefinite noise covariance matrix for 
all relevant wavenumbers (fc Q = 0, . . . , tt). Note that this 
region in general also depends on the relaxation rates A e 
and A e (we have chosen A e = A e = 1 in the plot), but 
it is independent of the temperature and density. The 
actually allowed combinations of simulation parameters 
are a subset of the white region in Fig. [21 as numerical 
stability provides further restrictions. In particular, a 
lower bound of the interface width leads, via eq. ([3]), to 
an upper bound on the sound speed in the liquid phase 
for each value of k (dashed curve in Fig. [5]) . 

We finally remark that the fluctuation temperature T, 
which is related to the velocity fluctuations by eq. ©, 
is bounded from above by the stability constraint of LB 
(known as low Mach-numbcr constraint in the ideal gas 
case), k B T/p = (u 2 ) < a 2 , or 



k B T < <j s p . 



(36) 
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FIG. 2: (Color online) Allowed region of simulation parame- 
ters. The white area corresponds to the combinations of ft and 
c s that ensure a positive-semidefinite noise covariance matrix 
H(k) for all k. The dashed red curve marks the parameter 
combinations that result in an interface width of 4 lattice 
units. Points above the curve correspond to smaller interface 
widths and thus can lead to potentially unstable simulations 
when interfaces are present. (All relaxation times are set to 
1.0.) 
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While this constraint also holds for simulations of ther- 
mal fluctuations in the ideal gas (24[, it becomes par- 
ticularly relevant in a two-phase system, as the above 
constraint must be fulfilled both in the liquid and gas 
phases with a uniform temperature throughout the sys- 
tem. Thus, in a two-phase system, the maximal attain- 
able temperature is limited by the vapor density. In 
particular, one has kgT = pz,(u| a ) = Pg( u g a ), hence 
(u 2 L Q ) = (pg/pl)(ug «)> showing that the effects of ther- 
mal fluctuations on the kinetics in the liquid phase reduce 
with increasing density ratio. This implies furthermore, 
that systems with high density ratios are effectively simu- 
lated at correspondingly larger length scales, where fluc- 
tuations are less pronounced. 



IV. IMPLEMENTATION 

We shall briefly describe here the practical steps re- 
quired to implement thermal noise in a simulation based 
on the modified-equilibrium model. Although the the- 
oretical analysis in the preceding section has been per- 
formed for an underlying D2Q9 lattice, the steps in the 
derivation of the FDT are nevertheless general and appli- 
cable to any DdQn lattice. The resulting noise covariance 
matrix 5 is expected to remain of similar form to (|34p . 

As mentioned in the preceding section, there exist in 
principle two options to model thermal fluctuations. The 
computationally easiest one is to employ spatially uncor- 
related noise, in which case the quantities £ a ( r ) hi the 
LBE (|23[) are spatially independent Gaussian random 
variables of covariance given by the zero wavenumber- 
limit 3(0) = lim fc ^o S(k) of the noise of eq. ([53]). This 
form of noise has the advantage of being readily applica- 
ble to inhomogeneous systems — however, at the expense 
of accepting equilibration errors at higher wavenumbers 
unless simulations are run in a rather small range of pa- 
rameters, as will be demonstrated in section [Vj If, on 
the other hand, only the correct behavior of the model 
at the largest length scales (i.e. the hydrodynamic limit) 
is of interest, the use of spatially uncorrelated noise is 
sufficient. Although it fulfills the FDT of fluctuating hy- 
drodynamics strictly only for k — > 0, one can expect it to 
give still acceptable results almost up to k ~ 1. 

The fluid momentum in a simulation subject to ther- 
mal noise should obey Gaussian statistics in real space 
with a variance given by (jajp) = PaksTSa^- Hence, 
as a first check in simulation whether uncorrelated noise 
is sufficient, one can track the momentum variance com- 
puted over the lattice in real-space and compare with the 
theoretically expected result. Since violated equiparti- 
tion at smaller scales will inevitably show up in the glob- 
ally computed equal-time variance, a deviation from the 
expected value indicates a length-scale dependent dissi- 
pation mechanism not captured by spatially uncorrelated 
noise. Consequently, the use of correlated noise would be 
required. More detailed assessments can be performed 
by computing wavcnumber-resolved variances (see sec- 



tion EJ. 

Due to the existence of cross-correlations between 
stress and ghost noises, indicated by a non-zero 349, the 
noise construction becomes slightly more involved com- 
pared to the ideal gas case, where the covariance ma- 
trix is diagonal. However, since the dynamics of ghost 
modes arc immaterial for the hydrodynamic limit of a 
LB model, one might argue that in this limit, one can 
drop the off-diagonal components of 5 and proceed us- 
ing "ideal gas-like" noise, with variance defined by the 
diagonal elements of 3. Simulations have indeed shown 
that this is a feasible option, leading to satisfactory equi- 
libration of all LB modes for small wavenumbers. In this 
work, however, we shall stick to the exact expression for 
3. In this case, noise variables £ a that have the required 
(non-diagonal) covariance 5(0) can be constructed inde- 
pendently on each lattice site by a standard method j55j , 
which consists of Cholesky-factorizing the noise matrix 
as 3(0) = LL T , with a left-triangular matrix L. The 
noise variables follow as £ a = Lv a , where the v a are a set 
of independent Gaussian random variables of unit vari- 
ance. Note that, since 5(0) is only positive semi-definite, 
the Cholesky factorization has to be performed effectively 
with the lower-right (n — d—1) x (re — d — 1) block matrix 
of 5. Alternatively, one can compute the Jordan decom- 
position of 5, that is 5 = Sdia.g(ei)S T , where S is the 
orthogonal matrix of the eigenvectors of 5 and the are 
the eigenvalues. The required matrix L is then given by 
L = S diag( v / e7). 

In cases where uncorrelated noise leads to strong vio- 
lations of equipartition or where the correct behavior of 
the model at smaller scales is important, one must use 
spatially correlated noise, i.e. implement the exact FDT 
([Ml) at each point in A:-spacc. The algorithmic construc- 
tion of the noise can be performed in Fourier space, as 
described, for example, in [5(|. Here, one Fourier trans- 
form back to the real lattice space is required per time- 
step. The noises £ a (k) at a point k are computed using 
the same procedure described above for k = 0. 



V. RESULTS 

A. Equilibration tests 

As a first benchmark test, we check whether the noise 
defined by eq. (j34f leads to the correct equilibration of all 
the LB modes in a simulation. For this purpose, we per- 
form simulations in a two-dimensional, periodic, homoge- 
neous one-phase system of size 128 x 128 l.u., using either 
the full (spatially correlated) noise or its k — (uncor- 
related) approximation. All parameters are chosen such 
that numerical stability is ensured also if a two-phase in- 
terface were present. This imposes, in particular, a lower 
bound of approximately 4 l.u. on the interface width, 
eq. We use a fluctuation temperature T = 10~ 7 (set- 
ting kg — 1) for all simulations, which lies well within 
the intrinsic stability constraint (|36[) of LB. The magni- 
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tilde of the resulting velocity and density fluctuations is 
then of the order of 10 ~ 3 . All relaxation times are set 
to a value of r = 1.0. For a linearized, homogeneous 
non-ideal fluid model, the essential input parameters are 
the square-gradient parameter k and the thermodynamic 
speed of sound c s (which is equal to the generalized speed 
of sound in the limit k 0). The actual values of the 
parameters p^, py and /3 used in the bulk free energy 
density ([2]) arc immaterial in the linearized case, but are 
also stated for reference. 

Results are analyzed by comparing the equal-time cor- 
relations (|<5m a (k)| 2 ) of a LB mode to the theoretical 
expectation as expressed by the correlation matrix G, 
eq. (|33[) . This comparison is most conveniently per- 
formed by computing the equilibration ratio, which is 
defined as (|<5m a (k)| 2 ) sim /G aa (k). In the plots below, 
this quantity is shown as an average over 400 simulation 
snapshots. 

In Fig. [3j the equilibration ratios obtained with the 
exact (spatially correlated) form of noise are shown for 
K = 0.08 and c s = 0.27 (with the corresponding param- 
eters in the bulk free energy being py = 0.5, ph = 1.0, 
(3 = 0.14). We see that the maximum error in the equi- 
libration of every mode stays always below 5%, even for 
the largest wavenumbers. Quantitatively very similar re- 
sults have been obtained for all tested combinations of k, 
c s and r. In all cases, the equilibration error at interme- 
diate and large wavenumbers is found to remain less than 
10% (occasionally 20%), while it is generally negligible for 
smaller wavenumbers where all modes appear perfectly 
equilibrated. These results suggest that the noise covari- 
ance matrix (|34[) . and hence the underlying equilibrium 
correlation matrix G (f3"3"|) . correctly describe the dissi- 
pation and the fluctuations in the modified-equilibrium 
model. 

In contrast, using spatially uncorrelated noise with the 
above choice of simulation parameters has been found to 
lead to large errors at higher wavenumbers. However, at 
least for certain choices of k and c s , satisfactory results at 
all wavenumbers can be achieved also with uncorrelated 
noise. This is demonstrated in Fig. [4j where equilibra- 
tion ratios obtained for k = 0.03 and c s = 0.15 are shown 
(with the corresponding parameters in the bulk free en- 
ergy © chosen as p v = 0.1, p L = 1.0, /3 = 0.015). Wc 
see that the maximum error in the equilibration of the 
density and momentum modes is always below 10%, even 
for the largest wavenumbers. The errors in the stress 
and ghost modes, however, are significantly larger, espe- 
cially for intermediate k. These deviations are found to 
diminish if spatially correlated noise is used. In the hy- 
drodynamic region, which is located between k = and 
~ 0.8 for the present choice of parameters (57[, all errors 
can be considered as negligible. This is found to be true 
for all tested parameter combinations, and holds even for 
those cases where no acceptable thcrmalization at larger 
wavenumbers can be obtained with uncorrelated noise. 



B. Capillary fluctuations 

In the presence of interfaces, the fluctuations in the 
bulk fluid can induce capillary (or interfacial height) fluc- 
tuations [58[ . The static spectrum of large- wavelength, 
small-am plit ude fluctuations of the interface height h is 
given by US [H 

(|Mk)| 2 ) = ^, (37) 
<j k 

where a is the surface tension, eq. (01 . In order to test 
whether this relation can be reproduced by the fluctu- 
ating non-ideal fluid model, we perform simulations of a 
liquid-vapor interface that belongs to an extended liquid 
stripe placed in a fully periodic, two-dimensional box of 
size L x x L y = 2048 x 400 l.u. The liquid stripe has 
a width of 180 l.u. and is aligned parallel to the £-axis. 
Simulation parameters are pl = 1-0, pv = 0.5, ft = 0.04, 
k = 0.03, T = 10~ 7 and r = 0.1. For the simulation 
of capillary fluctuations with the modified-equilibrium 
model, Galilean-invariance correction terms become im- 
portant; we use the expressions proposed by Js3|- The 
interface height h(x) is defined as the position of the 
point where the density equals (pl + Pv)/%- Spatially 
uncorrelated Langevin noise defined by 5(0), cq. (|34[) . is 
employed. In order to apply this noise to inhomogeneous 
situations, such as the present one, the local values of the 
density po and the speed of sound c s have to be taken into 
account in computation of the real-space noise variance. 

Fig. E] shows the spectrum of the interfacial height fluc- 
tuations obtained for the modified-equilibrium model. 
We see, that the simulation results are well described 
by the theoretical capillary structure factor (|37|) for 
wavenumbers k < 0.2 (6(|. We ascribe the deviations 
at larger wavenumbers to the presence of error terms in 
the hydrodynamic equations of the model and to the fact 
that the harmonic approximation, on which expression 
(|37|) is based, breaks down and curvature corrections as 
well as lattice effects become important. At the smallest 
wavenumbers, fluctuations are damped due to compress- 
ibility of the liquid stripe whole. 

VI. SUMMARY AND OUTLOOK 

We have presented in this paper a systematic study 
of thermal fluctuations in the LB method for non- 
ideal fluids. As a specific example, we considered the 
Langevin version of the modified-equilibrium model pro- 
posed by Swift et al., which approximates the stochas- 
tic Navier-Stokes equations for non-ideal fluids based on 
square-gradient free energy functionals. Thermal fluc- 
tuations are implemented by adding Gaussian random 
noise sources to the deterministic LBE, thereby promot- 
ing it to a LB-Langevin equation. In order to deter- 
mine the covariance matrix of the noise, a fluctuation- 
dissipation theorem has been derived, which requires as 
input an expression for the equal-time correlations of the 



11 



1.20 
1.15 
1.10 
1.05 
1.00 
0.95 
0.90 
0.85 
0.80 



(a) 



(d) 



0.0 0.5 



1.0 1.5 
k 



1.2 
1.0 
0.8 
0.6 
0.4 
0.2 
0.0 











On o' v 








♦ k=2.5 




■ k=1.5 




• k=0.5 





2.0 2.5 3.0 



1.2 
1.0 

ST 

5? 0.8 

^ 0.6 

ST 

Q. 0.4 
0.2 
0.0 




♦ k=2.5 \$$\ 
■ k=1.5 

• k=0.5 



(b) 




1.20 
1.15 
1.10 
1.05 
1.00 
0.95 
0.90 
0.85 
0.80 



(c) 



(f) 



- • k 




° h 








!"!x...xX v „xx 






X X X 



0.0 0.5 



1.0 



1.5 

k 



2.0 2.5 3.0 



1.20 
1.15 
1.10 
1.05 
1.00 
0.95 
0.90 
0.85 
0.80 



. □ q x 






o q y 






A £ 






B B 8 ° ° = 






« S H 


! — ffl — * 




A A 

A A A A A 




i a u 

A A 



0.5 



1.0 



1.5 
k 



2.0 2.5 



3.0 



FIG. 3: (Color online) Correlated noise in the modified-equilibrium model for k = 0.08, c s = 0.265, r = 1.0. Equilibration 
ratios of (a,b) the density, (c,d) the momentum, (e) the transport and (f) the ghost modes. In (b) and (d) the dependence of 
the equilibration ratio of the density and momentum on 8, when k = (cos 8, sin9)k, is shown for several magnitudes of k. In 
the remaining plots, each data point represents an average over all directions in fc-space at each magnitude |fe|. j x denotes the 
x-component, j\\ the longitudinal and jt the transverse component (with respect to k) of the momentum j. 



distribution function in a non-ideal LB fluid. Drawing 
on continuum kinetic theory of fluids, a general ansatz 
[cq. (|27p] for these correlations is provided. Application 
of this ansatz to the modified-equilibrium LB model re- 
quires non-trivial modifications to ensure that the noise 
covariance matrix is positive-semidefinite and obeys the 
FDT of fluctuating hydrodynamics at finite wavenum- 
bers. We have obtained a wavenumber-dependent form 
of thermal noise that was shown to lead to excellent ther- 
malization of all modes in a simulation at all wavenum- 
bers and all tested parameters. The necessity of spatially 
correlated noise for the modified-equilibrium model, al- 
though a priori unexpected, can be traced back to the 
fact that this model is based on a non-Maxwellian form 
of the equilibrium distribution function, in contrast to 
standard kinetic theory. 

Additionally, we have demonstrated that thermal fluc- 
tuations in the hydrodynamic regime can often already 
be satisfactory modeled using noise evaluated in the 
zero wavenumber- limit. This result is expected to hold 
strictly whenever non-ideal fluid interactions enter the 
hydrodynamic equations reversibly, in which case the 
random stress tensor is independent of wavenumber. 
The modified-equilibrium model only approximately ful- 
fills this requirement, as the non-ideal interactions — even 
though weakly — do contribute to the dissipative terms in 
this model. Since noise obtained in the zero wavenumber- 
limit is spatially uncorrelatcd by construction, it has 
the advantage of being easy to implement and read- 
ily applicable to inhomogeneous systems. In the small 



wavenumber-regime, good equilibration of all LB modes 
has been obtained for all simulation parameters. 

For a future study, it would be interesting to apply 
the presented approach also to force-based non-ideal fluid 
models. Preliminary results using the model of Lee and 
Fischer [5^] indicate that acceptable equilibration for low 
wavenumbers can be obtained u sing noise identical to the 
one employed for the ideal gas [57|. This finding agrees 
with the expectations based on the corresponding FDT 
in the hydrodynamic limit. 

Although we assumed that density fluctuations are de- 
scribed by a square-gradient free energy functional, the 
present theory puts in fact no constraints on the form 
of the structure factor. Hence, our approach should in 
principle also be applicable to models based on differ- 
ent thermodynamic approaches, such as [6(j. Along this 
direction, the extension of the present theory to multi- 
component approaches such as, for example, binary fluids 
[29l ] . or emulsions 61 might be particularly interesting. 
Finally, it might be interesting to investigate how the 
noise for a non-ideal LB fluid can be derived without re- 
sorting to results of continuum kinetic theory, but instead 
using the approach to the fluctuating LBE proposed by 
Diinweg et al. [25} . 



VII. ACKNOWLEDGMENTS 

We thank Kevin Stratford, P.T. Sumcsh and Alexan- 
der Wagner for useful discussions. M.G. thanks the 



12 



1.20 
1.15 

3 i.io 

in 

X 105 
1.00 

3 0.95 

^ 0.90 

^ 0.85 

0.80 



0.0 0.5 1.0 



(a) 



1.5 2.0 
k 



1.2 
1.0 

H 

§ 0.8 
C 0.6 

s 

■A 0.4 
0.2 
0.0 











— \ 




' ♦ k=2.5 \^«°'\ ' 
, ■ k=1.5 




.k=0.5 ^ 


J J_ 



2.5 3.0 



1.2 
1.0 

ST 

5? 0.8 

^ 0.6 

ST 

Q. 0.4 
0.2 
0.0 







~ 




' ♦ k=2.5 " 
, ■ k=1.5 
• k=0.5 


o r 



1.20 r 



1.15 * Jx 

H o 1.10 \ °Jj 
~- 1.05 



1.00 
0.95 
0.90 
0.85 
0.80- 



-i-S- 



•» * s » jj s ° 



« if * 



0.0 0.5 



1.0 



(b) 



(c) 



1.5 

k 



2.0 2.5 



3.0 



(d) 




1.3 



8 1.2 



1.1 



1.0 
0.9 



□ q x 
o q. 



0.5 



1.0 



(f) 



1.5 
k 



2.0 



2.5 



3.0 



FIG. 4: (Color online) Uncorrelated noise in the modified-equilibrium model for k = 0.03, c s = 0.15, r = 1.0. Equilibration 
ratios of (a,b) the density, (c,d) the momentum, (e) the transport and (f) the ghost modes. In (b) and (d) the dependence of 
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Appendix A: Kinetic theory of non-ideal fluids 

We present here a brief account on the kinetic theory 
of fluctuations in a non-ideal fluid and show how this the- 
ory can be applied to the LB method and the particular 
model considered in section UlI CI 



1. Continuum theory 

In a kinetic description of a non-ideal fluid, the knowl- 
edge of the density and momentum correlations [eqs. ([?]) 
and ©] alone is not sufficient. Instead, one must spec- 
ify also the equal-time correlations of the fluctuations in 
the one-particle distribution function /(r,c), where c is 
the molecular velocity. The corresponding expression can 
be motivated from statistical mechanics within Klimon- 
tovich's approach to kinetic theory fl9L |62|. 

This approach is based on defining a phase-space den- 
sity as 

F(r, c) = fJ.y^5(r - r q )6(c - a) , 
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where p, r and c refer to the mass, positions and veloci- 
ties of the fluid particles, respectively. Since in our case 
the relevant quantity is the mass density instead of the 
number density, we have defined the phase density with 
an additional factor of mass. This ensures that the one- 
particle distribution functions derived below are defined 
in terms of mass density, in agreement with the situation 
in LB. Reduced particle distribution functions f n can be 
defined by computing moments of F with respect to the 
full TV-particle distribution function / . The first two 
reduced distribution functions are given by [62j : 



2. LB theory 



(F(r ! c))=/ 1 (r ! c) 
( J F(r,c)F(r',c / )>=/ 2 (r ! c,r',c / ) 

+ /!<S(r-r')<5(c-c')/i(r,c) 



(Al) 
(A2) 



We consider now fluctuations SF = F — (F) around a 
global equilibrium state with density po and zero flow 
velocity, i.e. (F) = /i(c) is assumed to be a global 
Maxwcllian distribution. The fluctuation SF can then 
be interpreted as a fluctuation Sfi in the one-particle 
distribution function f\ over a uniform reference state 
described by f\. In a Langevin description, we thus pro- 
mote /i to be an instantaneously fluctuating quantity 
given by /i(r,c) = /i(c) + <5/i(r,c). Specializing to a 
translationally invariant system, eq. (|A2[) allows to de- 
termine the equal-time correlations of the fluctuations in 
the one-particle distribution function as 

(*/i(r,c)*/i(r , ,c')) = (SF(r, C )SF(r',c')) 

= (F(r,c)F(r',c'))-/i(c)/ 1 (c') 

= /i(c)/i(c')[.9(r-r')-l] 
+ //S(r-r')5(c-c')/i(c), 

(A3) 

where in the last step, we introduced the pair correlation 
function g by / 2 (r - r', c, c') = / 1 (c)/ 1 (c').g(r - r'). The 
structure factor, defined in terms of the relative fluctu- 
ations, S(r) = (Sp(r)Sp) [cf. Q], is related to the pair 
correlation function by [9|, l37| 

S(r) = pl(g(v)-l)+pp Q S(r). 



Transforming relation (IA3[) to Fourier space and drop- 
ping the index 1, we finally obtain the desired relation 

(J/(k,c)5/(k',c')) = [f(c)f(c')[S(k)/ PQ -p}/p 

+ pf(c)S{c - c')l S(k + k') . (A4) 



The first term on the right hand side of (|A4[) describes 
spatial correlations due to the non-ideal character of the 
fluid. For the ideal gas, S(k) = ppo, thus only the last 
term remains, and relation (|A4[) becomes identical to the 
expression used in the Boltzmann-Langevin theory for a 
dilute gas [l|| [l?], [H| . As can be easily checked, relation 
(|A4[) contains already expression ([7]) for the structure 
factor and expression ([5]) for the momentum correlations, 
since, by definition, p = J f(c)d d c and j — J f(c)c d d c. 



The LB analog of eq. (|A4p for the equal-time corre- 
lations of fluctuations in the distribution function of a 
non-ideal fluid, can be written as [eq. (|2"T|)] 

(SMk)Sf 3 (k')) = [fifjlSW/po-A/po+tfiSi^V&k,-* , 

. (A5) 

where fi = /j Cq (poi u = 0) is the distribution function of 
a quiescent reference state. In this context, the parame- 
ter p can be interpreted as the mass of a fictitious fluid 
particle. Its value is a priori unknown and thus has to 
be determined such that the correct momentum variance 
as required by statistical mechanics, eq. ([5]), is obtained. 
The factor V, representing the system volume, arises by 
dimensional consistency and will be neglected henceforth. 
The correlation matrix G is obtained from (|A5[) through 
a basis transformation, 



G ab (k) =T a ir bj (Sfi(k)Sf j (-k)) 

= fh a fh b [S(k)/p - p]/p + pT al T bi fi 



(A6) 



We shall briefly demonstrate how the unknown param- 
eter p can be determined for the example of an ordinary 
ideal gas LB model. In this case, we have fi = poWi, 
Y^i c ifi = and hence 



(jajp) = C la C k j3(5f l 8f k ) 



*Cif3pfi = PP0<J^S a j3 



where additionally the orthogonality of the basis vectors 
and relation (|2"0|) have been invoked. In order to obtain 
the desired expression {jajp} = PokBTS a p, we see that 
the mass parameter must be chosen as p = ksT/a^ 
in agreement with the ideal gas equation of state. An 
analogous calculation of the density correlator shows that 
p is in fact the structure factor of the ideal LB gas divided 
by Po, P = SidXB/Po- 

A strict application of eq. (|A6[) to compute the 
equilibrium correlation matrix G for the modified- 
equilibrium model requires to use the equilibrium mo- 
ments (Table [TTJ evaluated in a quiescent state, fh a = 
m^(p 0> u=0) = {po^0,d(0)po,0,0,0,0,-d(0)p }, 
where d(0) = lim^o d(k). The requirement that G\\ = 
S(k) and G 22 = G 33 = Pok B T fixes p = k B T/c 2 s . 
However, it can be shown that if the equilibrium cor- 
relation matrix G constructed in this way is used, the 
noise following from eq. (|26|) violates the FDT fluctu- 
ating hydrodynamics at any finite wavenumber. This 
can be understood from the fact that the bulk viscosity 
of the modified-equilibrium model is given by eq. (|35p . 

C(k) = p a1 [r b — |) 2 — , and — in particular — is 

dependent on wavenumber. As shown in appendix[Cl the 
FDT of fluctuating hydrodynamics requires the same fac- 
tor 2 — c^(k)/tTg to appear in variance of the LB noise 
pertaining to the bulk stress mode, (£f). The above G, 
however, would lead to an expression for (£f ) that is cor- 
rect only at k = 0. Moreover, a closer investigation of 
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the so obtained noise matrix reveals the presence of a 
negative eigenvalue for finite wavenumbers — a fact that 
invalidates its meaning as a covariancc matrix and in- 
hibits its use in simulation. 

All the above mentioned problems can, however, be 
successfully solved by replacing all occurrences of c 2 in 
G(k) by the full fc-dcpcndent speed of sound c^(k) = 
pok B T/ S(k). This is tantamount to re-introducing the 
/c-dependent terms c?(k) in m a , and then using the so 
defined fc-dependent reference state m a (k) in the com- 
putation of G according to (| A6|) . Agreement with the 
statistical mechanical expression for the momentum cor- 
relation function now requires to introduce a fc-dependent 
"mass" parameter p as 

M (k) = fc B r/c2(k) = 5(k)/ Po . 



Inserting this expression for p(k) into (| A5|) then finally 
leads to relation (|32|) . 

The above derivation shows that, while equation (|A5[) 
is a useful starting point to arrive at a valid expression 
for the correlations of the distribution function, it may 
have to be modified when applied to a specific non-ideal 
LB fluid. In the case of the model of Swift et al. we 
considered above, the modifications can be traced back to 
the fact that the local equilibrium distribution is defined 
in a non-standard way compared to continuum kinetic 
theory, where one assumes a Maxwcllian form. Hence, 
we expect relation (|A5[) in its original form to be more 
appropriate to non-ideal fluid LB models that employ 
the usual ideal gas equilibrium distribution. This will be 
investigated in future works. 



Appendix B: Lattice Fourier transforms 



text. In general, deviations between the continuum and 
discrete derivative operators become significant only at 
intermediate and high wavenumbers. There, isotropy is 
typically lost, i.e., the lattice Laplace operator depends 
on the direction in fc-space. 

Appendix C: Fluctuating stress tensor 

In the hydrodynamic limit, it is possible to directly 
compute the noise strength of the stress modes required 
by the FDT of fluctuating hydrodynamics, cqs. (fT6|) and 
(|16l) . In Fourier space, the random stress tensor correla- 
tions can be expressed as (ioj 



(i? a/ 3(k,t)i? 7 5(k',t')) =2k B T rj[8 ai 5f}8 + S a s 



- -5 a p5 lS j + C5 a p5 1 s <5k,-k"5t,i', (CI) 

where the viscosities are in principle allowed to depend on 
the magnitude of the wavenumber k. On the other hand, 
a Chapman-Enskog analysis performed on the fluctuating 
LBE (f2"3"|) reveals that, for the presently chosen D2Q9 
basis set (see Table IJ), the random stress tensor is related 
to the noise variables by (cf. [25j]) 



?4 



is 



R, 



a/3 



\ b VNE X, 
He 



XsVNe 
Ji £5 



(C2) 



/NE 



XsVNE, 



This result is independent of the particular non-ideal 
fluid model and only depends on the underlying lattice. 
After rearranging (|C2[) to obtain the £ a in terms of R a p , 
we use (|Clj) to evaluate the variances and finally plug 
in the known expressions of the shear and bulk viscosity. 
We find 



Working with Fourier transformations on a lattice re- 
quires to take into account the correct equivalents of the 
discretized derivative operators. In the present case, a 
discretized Laplace operator of the form 



[p(r + a) + p(r - a) - 2p(r)] /a 3 . 



is used, where i runs over all eight non-zero directions of 
the D2Q9 lattice. The Fourier transformed discretized 
Laplace operator follows as 



'4 

— (cosk x + COS ky) ■ 



2 

-(cos k. 



This expression has to be used in place of — k 2 in the LB 
analogs of the structure factor and the speed of sound. 
Note that, in order not to clutter- up notation, we prefer 
to state the continuum expressions throughout the main 



£4 = 3A b Tri? 
(Cf) = 9\ 2 b ■ 8k B TC 

= -36k B Tpa 2 (2X b + X 2 b )h, 



(C3) 



where h = 1 for an ideal gas-like model and 
h = 2 — c 2 (k)/a 2 for the modified equilibrium model due 
to the modified bulk viscosity ([55)1 . 



£5 — ^s{Rxx — Ryy) 

(g) = Xl ■ 8fc B Tr ? 

= -Ak B Tpa 2 s {2\ s + \ 2 s ) 



^6 — X S (R X y + Ry X )/2 

(4 2 ) = X 2 S ■ 2k B T V 
= -k B Tpa 2 {2X s 



xl). 



[1] S. Succi, The Lattice Boltzmann Equation for Fluid Dy- [2] M. Rauscher and S. Dietrich, Annu. Rev. Mater. Res. 38, 
namics and Beyond (OUP, Oxford, 2001). 



15 



143 (2008). 

D. Bonn, J. Eggers, J. Indekeu, J. Meunier, and E. Rolley, [36 
Rev. Mod. Phys. 81, 739 (2009). [37 
P. C. Hohenberg and B. I. Halperin, Rev. Mod. Phys. 49, 

435 (1977). [38 
A. J. Bray, Adv. Phys. 43, 375 (1994). 
R. F. Fox and G. E. Uhlenbeck, Phys. Fluids 13, 1893 [39 
(1970). 

E. H. Hauge and A. Martin-Loef, J. Stat. Phys. 7, 259 [40 
(1973). 

A. J. C. Ladd, Phys. Rev. Lett. 70, 1339 (1993). [41 
P. M. Chaikin and T. C. Lubensky, Principles of Con- 
densed Matter Physics (Cambridge, 1995). [42 
J. M. O. de Zarate and J. V. Sengers, Hydrodynamic 
Fluctuations in Fluids and Fluid Mixtures (Elsevier, [43 
2006). [44 
K. Kadau, C. Rosenblatt, J. L. Barber, T. C. Germann, 

Z. Huang, P. Carles, and B. J. Alder, Proc. Nat. Acad. [45 
Sci. 104, 7741 (2007). 
M. Moseler and U. Landman, Science 289, 1165 (2000). [46 

B. Davidovitch, E. Moro, and H. A. Stone, Phys. Rev. 

Lett. 95, 244505 (2005). [47 

R. Fetzer, M. Rauscher, R. Seemann, K. Jacobs, and [48 

K. Mecke, Phys. Rev. Lett. 99, 114503 (2007). 

B. B. Kadomtsev, Sov. Phys. JETP 5, 771 (1957). [49 

M. Bixon and R. Zwanzig, Phys. Rev. 187, 267 (1969). 

R. F. Fox and G. E. Uhlenbeck, Phys. Fluids 13, 2881 [50 

(1970). 

L. D. Landau and E. M. Lifshitz, Fluid Mechanics (Perg- [51 
amon, 1959). 

Y. L. Klimontovich, Sov. Phys.-Usp. 16, 512 (1973). [52 
B. Kim and G. F. Mazenko, J. Stat. Phys. 64, 631 (1991). 
A. J. C. Ladd, J. Fluid Mech. 271, 285 (1994). [53 
A. J. C. Ladd, J. Fluid Mech. 271, 311 (1994). 

A. J. C. Ladd and R. Verberg, J. Stat. Phys. 104, 1191 [54 
(2001). 

R. Adhikari, K. Stratford, M. E. Cates, and A. J. Wagner, [55 
Europhys. Lett. 71, 473 (2005). 

B. Duenweg, U. D. Schiller, and A. J. C. Ladd, Phys. [56 
Rev. E 76, 036704 (2007). 
J. W. Dufty and M. H. Ernst, in Pattern Formation [57 
and Lattice Gas Automata, edited by A. Lawniczak and [58 
R. Kapral (Am. Math. Soc, 1993). [59 
J. B. Bell, A. L. Garcia, and S. A. Williams, Phys. Rev. [60 
E 76, 016708 (2007). [61 
M. R. Swift, W. R. Osborn, and J. M. Yeomans, Phys. 

Rev. Lett. 75, 830 (1995). [62 
M. R. Swift, E. Orlandini, W. R. Osborn, and J. M. [63 
Yeomans, Phys. Rev. E 54, 5041 (1996). 

J. S. Rowlinson and B. Widom, Molecular Theory of Cap- [64 
illarity (Dover Publications, 1982). 

D. Jamet, O. Lebaigue, N. Coutris, and J. M. Delhaye, [65 
J. Comp. Phys. 169, 624 (2001). 
D. M. Anderson, G. B. McFadden, and A. A. Wheeler, [66 
Annu. Rev. Fluid. Mech. 30, 139 (1998). 
A. J. M. Yang, P. D. F. Ill, and J. H. Gibbs, J. Chem. 
Phys. 64, 3732 (1976). 
R. Evans, Adv. Phys. 28, 143 (1979). 
H. Goldstein, Classical Mechanics (Addison- Wesley, 



1980). 

Q. Zou and X. He, Phys. Rev. E 59, 1253 (1999). 
J.-P. Hansen and I. R. McDonald, Theory of Simple Liq- 
uids (Academic Press, 2006), 3rd ed. 
L. D. Landau and E. M. Lifshitz, Statistical Physics. Part 
I (Pergamon, 1980). 

J. P. Boon and S. Yip, Molecular Hydrodynamics 
(McGraw-Hill, 1980). 

L. E. Reichl, A Modern Course in Statistical Physics 
(Wiley, 1998), 2nd ed. 

A.-M. S. Tremblay, M. Arai, and E. D. Siggia, Phys. Rev. 
A 23, 1451 (1981). 

H. Risken, The Fokker-Planck Equation (Springer, 1989), 
2nd ed. 

J. Keizer, Phys. Fluids 21, 198 (1978). 

F. J. Higuera and J. Jimenez, Europhys. Lett. 9, 663 

(1989). 

F. J. Higuera, S. Succi, and R. Benzi, Europhys. Lett. 9, 
345 (1989). 

R. R. Nourgaliev, T. N. Dinh, T. G. Theofanous, and 
D. Joseph, Int. J. Multiphase Flow 29, 117 (2003). 
D. Raabe, Model. Simul. Mater. Sci. Eng. 12, R13 (2004). 
R. Benzi, S. Succi, and M. Vergassola, Phys. Rep. 222, 
145 (1992). 

C. K. Aidun and J. R. Clausen, Annu. Rev. Fluid Mech. 
42, 439 (2010). 

D. d'Humieres, Prog. Astronaut. Aeronaut. Ser. 159, 459 
(1992). 

R. Adhikari and S. Succi, Phys. Rev. E 78, 066701 
(2008). 

P. Lallemand and L.-S. Luo, Phys. Rev. E 61, 6546 
(2000). 

C. M. Pooley and K. Furtado, Phys. Rev. E 77, 046702 
(2008). 

D. J. Holdych, D. Rovas, J. G. Georgiadis, and R. O. 
Buckius, Int. J. Mod. Phys. C 9, 1393 (1998). 

W. H. Press et al., Numerical Recipes (Cambrigde Uni- 
versity Press, 2007), 3rd ed. 

J. Garcia-Ojalvo and J. M. Sancho, Noise in Spatially 
Extended Systems (Springer, 1999), 1st ed. 
M. Gross et al., unpublished (2010). 
M. Grant and R. C. Desai, Phys. Rev. A 27, 2577 (1983). 
T. Lee and P. F. Fischer, Phys. Rev. E 74, 046709 (2006). 
X. Shan and H. Chen, Phys. Rev. E 47, 1815 (1993). 
R. Benzi, S. Chibarro, and S. Succi, Phys. Rev. Lett. 
102, 026002 (2009). 

R. L. Liboff, Kinetic Theory (Springer, 2003), 3rd ed. 
L. D. Landau and E. M. Lifshitz, Physical Kinetics (Perg- 
amon, 1981). 

We neglect here any dimensional factors due to delta 
functions 5(0) evaluated at zero. 

k 2 has to be understood here as the negative of the 
Fourier-transformed Laplace operator, see appendix [Bl 
Simulations using ideal gas- type noise for the force-based 
model of Lee and Fischer [5(3] gave similar results, agree- 
ing with the capillary structure factor up to a slightly 
larger value of k ~ 0.5. 



